#!/bin/perl
# Author Edgar Medina
use Bio::SeqIO;
use Bio::DB::Fasta;
use Getopt::Long;
use strict;
use warnings;

#script for getting the Ids from HMMER
#usage: perl do_indexngetseq2.pl --dbname DIRECTORYOFSEQUENCES HMMEROUT.TAB > RETRIEVEDSEQUENCES
my $dbname;
my $cutoff = 1e-10;
GetOptions(
	'c|cutoff:s' => \$cutoff,
	'dbname:s'   => \$dbname,
	 );

my %unique;
while (<>) {
	next if /^\#/;
	chomp;
	my ($target,$targetacc,$query,$queryacc,$evalue,$score) = split(/\s+/,$_);
	next unless $evalue <= $cutoff;
	$unique{$target}++;
}

#SCRIPT FOR INDEXING A LOCAL DATABASE PRIOR TO FETCHING SEQUENCES!!!!
# if you provided a database...(dbname)
# Index it
if ( $dbname ) {
my $db = Bio::DB::Fasta->new($dbname);
my $out = Bio::SeqIO->new(-fh => \*STDOUT, -format => 'fasta'); # our SeqIO fasta writer object

 # cleaned ID, let us fetch it from the database
for my $id ( keys %unique ) {
 my $seq = $db->get_Seq_by_id($id);
 if( $seq ) {
   my $temp_id = $seq->id;
   # could fix the ID to remove any prefixes if you want
   # maybe more checks on the ID and cleaning it up next
   if( $temp_id =~ /^jgi\|/) {
	my $tseq = Bio::Seq->new(-id => $seq->id,
                            -desc => $seq->description,
                            -seq => $seq->seq);
	my @row = split(/\|/,$temp_id); 
   	 $seq = $tseq;
        $temp_id = $row[1] . '_' . $row[-1]; # -1 means the LAST one 
        $seq->id($temp_id);
   }
   $out->write_seq($seq); 
 } else {
   warn("Cannot find sequence '$seq' in the db, did you parse the ID properly?\n");
 }
}
}
# put subroutines here

